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Abstract 

The increasing size of cosniological simulations has led to the need for new visualization techniques. 
We focus on Smoothed Particle Hydrodynamical (SPH) simulations run with the GADGET code 
and describe methods for visually accessing the entire simulation at full resolution. The simulation 
snapshots are rastered and processed on supercomputers into images that axe ready to be accessed 
through a web interface (GigaPan). This allows any scientist with a web-browser to interactively 
explore simulation datasets in both in spatial and temporal dimensions, datasets which in their native 
format can be hundreds of terabytes in size or more. We present two examples, the first a static 
terapixel image of the MassiveBlack simulation, a P-GADGET SPH simulation with 65 billion particles, 
and the second an interactively zoomable animation of a different simulation with more than one 
thousand frames, each a gigapixel in size. Both are available for public access through the GigaPan 
web interface. We also make our imaging software publicly available. 
Subject headings: Cosmology: observations - large-scale structure of Universe 



1. INTRODUCTION 

In the 40 years that iV-body simulations have been 
used in Cosmology research, visualization has been the 
most indispensable tool. Physical processes have often 
been identified first and studied via images of simula- 
tions. A few examples are: formation of filam entary 



structures in the large-scale distribution of matter ( Jenk 
ins et al IfrgM; Colberg et aT'SOOS'; Springel et al. 2UU5 
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Dolag et al. 2006), growth of feedback bubbles a rouni 
quasars (jsijacki et al.||2007| |Di Matteo et alj l2005|); cold 



fiows of ga s forming galaxies (Dekel et al. 2009 Keres 



et al. 



20051, and the evoluti on of ionization fronts dur 



the re-ionization epoch (Shin et al. 
|2007|) . The size of currem 



2008 
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Zahn et 
|al.||2007|. The size of current and upcoming Jr'eta-scale 
simulation datasets can make such visual exploration to 
discover new physics technically challenging. Here we 
present techniques that can be used to display images 
at full resolution of datasets of hundreds of billions of 
particles in size. 

Several implementations of visualization software for 
cosm ological simulations already exist. IRFIT (Gnedin] 
20111 is a general purpose visualization suite tnat can 
deal with mesh based scalar, vector and tensor data, as 
we ll as p article based datasets as points. YT (Turk et 
al. 20111 is an analysis toolkit for mesh based simula- 



tions that also supports imaging. SPLASH tPricejj2007j) 
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is a visualization suite specialized for simulations that 
use smoothed particle hydrodynamics(SPH) techniques. 
Aside from the CPU based approaches mentioned above, 
Szalay et al. ( 2008 1 implemented a GPU based interac- 

tive visualization tool for SPH simulations^ 

The Millennium I & II simul ations (^Springel et al. 



2005 Boylan-Kolchin et al.|2009 l have been used to test 
an interactive scalable rendering system developed by 
Fraedrich et al.| ( |2009[ ). Both SPLASH and the Mil- 
lenmum visualizer support high quality visualization of 
SPH data sets, while IRFIT treats SPH data as discrete 
points. 

Continuing improvements in computing technology 
and algorithms are allowing SPH cosmological simula- 
tions to be run with ever increasing numbers of parti- 
cles. Runs are now possible on scales which allow rare 
objects, such as quasars to form in a large simulation vol- 
ume with uniform high resoluti o n (see Section 2.1; and 
Di Matteo et al.| ([2011, in prep [); [Degraf et all ( |2011, in 
prep ~ Khandai et al. 1 2011, in prep \J. Being able to 
scan through a vast volume and seek out the tiny regions 
of space where most of the activity is occurring, while 
still keeping the large-scale structure in context neces- 
sitates special visualization capabilities. These should 
be able to show the largest scale information but at the 
same be interactively zoomable. However, as the size of 
the datasets quickly exceeds the capability of moderately 
sized in-house computer clusters, it becomes difficult to 
perform any interactive visualizations. For example, a 
single snapshot of the MassiveBlack simulation (Section 
2.1) consists of 8192 files and is over 3 TB in size. 

Even when a required large scale high resolution image 
has been rendered, actually exploring the data requires 
special tools. The GigaPan collaboration [^ has essen- 
tially solved this problem in the context of viewing large 
images, with the GigaPan viewer enabling anyone con- 
nected to the Internet to zoom into and explore in real 
time images which would take hours to transfer in total- 

* http://www.gigapan.org 



ity. The viewing technology has been primarily used to 
access large photographic panoramas, but is easily ap- 
plicable to simulated datasets. A recent enhancement to 
deal with the time dimension, in the form of gigapixel 
frame interactive movies (GigaPan Time MachincP]) will 
turn out to give particularly novel and exciting results 
when applied to simulation visualization. 

In this work we combine an off-line imaging technique 
together with GigaPan technology to implement an in- 
teractively accessible visual probe of large cosmological 
simulations. While GigaPan is an independent project 
(uploading and access to the GigaPan website is publicly 
available), we release our toolkit for the off-line visual- 
ization as Gae psp^ a software pack a ,ge aimed speci fically 
at GADGET |Spmigel et~al]|2001[ |Springel|[2005l ) SPH 
simulations. 

The layout of our paper is as follows. In Section 2 
we give a brief overview of the physical processes mod- 
eled in GADGET, as well as describing two P-GADGET 
simulations which we have visualized. In Section 3 we 
give details of the spatial domain remapping we employ 
to convert cubical simulation volumes into image slices. 
In Section 4, we describe the process of rasterizing an 
SPH density field, and in Section 5 the image render- 
ing and layer compositing. In Section 6 we address the 
parallelism of our techniques and give measures of per- 
formance. In Section 7 we briefly describe the GigaPan 
and GigaPan Time Machine viewers and present exam- 
ples screenshots from two visualizations (which are both 
accessible on the GigaPan websites). 

2. SIMULATION 



The MassiveBlack simulation is the state-of-art SPH 
simul ation of a AC DM universe ( Di Matteo et al.|2011, in 
prepl) . P-GADGET was used to evolve 2 x 3200'' particles 



Ada ptive Mesh Refinement (AMR, e.g., [Bryan fc Nor 
"^ H 



19971)) and Smo othed Particle Hydrodynamics 
Monaghan ( |1992[ )) are the two most used schemes 
rymg out cosmological simulations. In this work 



man 

for carrymg 

we focus on the visualization of the baryonic matter in 

SPH simulations run with P-GADGET ( Springel 2005). 

GADGET is an SPH implementation, and P-GXdGET 
is a version which has been developed specifically for 
petascale computing resources. It simultaneously fol- 
lows the self-gravitation of a collision-less N-body sys- 
tem (dark matter) and gas dynamics (baryonic matter), 
as well as the formation of stars and super-massive black 
holes. Dark matter particles and gas particle positions 
and initial characteristics are set up in a comoving cube, 
and black hole and star particles are created according to 
sub-grid mo deling ( Hernquist & Springel 2003 , Di Mat- 
teo et al. 2008| ) Gas particles carry hydrodynamical prop- 
erties, such as temperature, star formation rate, and neu- 
tral fraction. 

Although our attention in this paper is limited to imag- 
ing properties the of gas, stars and black holes in GAD- 
GET simulations, similar techniques could be used to vi- 
sualize the dark matter content. Also, the software we 
provide should be easily adaptable to t he data formats 
of ot her SPH codes (e.g. GASOLINE, ( [Wadsley et~ar 
2004| ) 



2.1. MassiveBlack 

® http://tiineniachine.gigapan.org 

^^ http://web.phys.cnm.edu/~yfengl/gaepsi 



in a volume of side length 533h^^Mpc with a gravita- 
tional force resolution of 5h~^Kpc. One snapshot of the 
simulation occupies 3 tera- bytes of disk space, and the 
simulation has been run so far to redshift z — 4.75, cre- 
ating a dataset of order 120 TB. The fine resolution and 
large volume of the simulation permits one to usefully 
create extremely large images. The simulation was run 
on the high performance computing facility, Kraken, at 
the National Institute for Computational Sciences in full 
capability mode with 98,000 CPUs. 

2.2. E5 

To make a smooth animation of the evolution of 
the universe typically requires hundreds frames directly 
taken as snapshots of the simulation. The scale of the 
MassiveBlack run is too large for this purpose, so we 
ran a much smaller simulation (E5) with 2 x 336^ par- 
ticles in a 50h~^Mpc comoving box. The model was 
again a AC DM cosmology, and one snapshot was output 
per 10 million years, resulting in 1367 snapshots. This 
simulation ran on 256 cores of the Warp cluster in the 
Mc Williams Center for Cosmology at CMU. 

3. SPATIAL DOMAIN REMAPPING 

Spatial domain remapping can be used to transform 
the periodic cubic domain of a cosmological simulation 
to a patch whose shape is similar to the domain of a sky 
survey, while making sure that the local structures in the 
simulation are preserved ( Carlson fc White|2010 Hilbert 
et al. 20091. Another application is making a thin slice 
that includes the entire volume of the simulation. Our 
example will focus on the latter case. 

A GADGET cosmological simulation is usually defined 
in the periodic domain of a cube. As a result, if we let 
f{X = (x,y,z)) be any position dependent property of 
the simulation, then 

f{X)^f{{x + iiL,y + yL,z + <jL)), 

where /i, v^ a are integers. The structure corresponds to 
a simple cubic lattice with lattice constant a — L, the 
simulation box side-length. A bijective mapping from 
the cubic unit cell to a remapped domain corresponds 
to a choice of the primitive cell. Figure [T| illustrates the 
situation in 2 dimensions. 

Whilst the original remapping algorithm by [Carlson fc| 
White (2010 1 results in the correct transformations being 
applied, it has two drawbacks: (i) the orthogonalization 
is invoked explicitly and (ii) the hit-testing for calcula- 
tion of the shifting (see below) is against non-aligned 
cuboids. The second problem especially undermines the 
performance of the program. In this work we present a 
faster algorithm based on similar ideas, but which fea- 
tures a QR decomposition (which is widely available as a 
library routine), and hit-testing against an A ABB (Axis 
Aligned Bounding Box). 

First, the transformation of the primitive cell is given 
by a uni-modular integer matrix, 

/Mil Mi2 Mi3\ 
M = M21 M22 M23 , 
VAfai M32 M33J 







Figure 1. Transformation of the Primitive Cell. The cubical unit 
cell is shown using solid lines. The new primitive cell, generated by 
[ j^ 2 ] '^ shown with dash-dotted lines. The transformed domain is 
shown in gray. 

where Mij are integers and the determinant of the matrix 
|M| = 1. It is straight-for ward to obtain such matrices 
via enumeration. ( [Carlson fc White 20101 The QR de- 
composition of M IS 

M = QR, 

where Q is an orthonormal matrix and R is an upper- 
right triangular matrix. It is immediately apparent that 
(i) application of Q yields rotation of the basis from the 
simulation domain to the transformed domain, the col- 
umn vectors in matrix Q-^ being the lattice vectors in 
the transformed domain; (ii) the diagonal elements of R 
are the dimensions of the remapped domain. For imag- 
ing it is desired that the thickness along the line of sight 
is significantly shorter than the extension in the other 
dimensions, thus we require < i?33 ^ |i?22| < |^ii|- 
Note that if a domain that is much longer in the line of 
sight direction is desired, for example to calculate long 
range correlations or to make a sky map of a whole sim- 
ulation in projection, the choice should be < |i?.33| < 

\R22\ « |i?ll|. 

Next, for each sample position X, we solve the in- 
definite equation of integer cell number triads / = 

X = Q'^X + aQ^I, 

where a is the box size, X is the transformed sample 
position satisfying X £ [0,i?ii) x [0, i?22) x [0, -R33). In 
practice, the domain of X is enlarged by a small number 
e to address numerical errors. Multiplying by Q on the 
left and re-organizing the terms, we find 



/ = 



QX 



X 

a 



Notice that QX is the transformed sample position ex- 
pressed in the original coordinate system, and is bounded 
by its AABB box. If we let (QX/a), G [Bi, T,], where Bi 
and Ti are integers, and notice — € [0, 1), the resulting 
bounds of / are given by 

l^e[B„T,]. 

We then enumerate the range to find X. 

When the remapping method is applied to the SPH 
particle positions, the transformations of the particles 
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Figure 2. Boundary effects and Smoothed Particles. Four images 
of a particle intersecting the boundary are shown. The top-right 
image is contained in the transformed domain, but the other three 
are not. The contribution of the two bottom images is lost. By 
requiring the size of the transformed domain to be much larger than 
typical SPH smoothing lengths, most particles do not intersect a 
boundary of the domain and the error is contained near the edges. 



Table 1 

Transformations 



Simulation 



MassiveBlack 



E5 



Matrix M 



[5 6 21 




[3 101 


3-7 3 




2 10 


L2 7 J 




Lo 2 1 J 



i?ii(h-iMpc) 



3300 



180 



Rii Pixels(Kilo) 



810 



36.57 



i?22(h-'Mpc) 



5800 



101 



i?22 Pixels(Kilo) 



1440 



20.48 



i?33(h-^Mpc) 



6.7 



Resolution(h ^Kpc) 



4.2 



4.9 



that are close to the edges give inexact results. The 
situation is similar to the boundary error in the origi- 
nal domain when the periodic boundary condition is not 
properly considered. Figure [2] illustrates the situation 
by showing all images of the particles that contribute 
to the imaging domain. We note that for the purpose 
of imaging, by choosing a i?33(the thickness in the thin- 
ner dimension) much larger than the typical smoothing 
length of the SPH particle, the errors are largely con- 
strained to lie near the edge. These issues are part of 
general complications related to the use of a simulation 
slice for visualization. For example in an animation of 
the distribution of matter in a slice it is possible for ob- 
jects to appear and disappear in the middle of the slice as 
they pass through it. These limitations should be borne 
in mind, and we leave 3D visualization techniques for 
future work. 

The transformations used for the MassiveBlack and E5 
simulations are listed in Table [TJ 

4. RASTERIZATION 

In a simulation, many field variables are of interest in 
visualization. 

• scalar fields : density p, temperature T, neutral 
fraction xm, star formation rate 0; 

• vector fields: velocity, gravitational force. 



In an SPH simulation, a field variable as a function of 
spatial position is given by the interpolation of the par- 
ticle properties. Rasterization converts the interpolated 
continuous field into raster pixels on a uniform grid. The 
kernel function of a particle at position y with smoothing 
length h is defined as 



Wiy,h) = 



r/i3 



1-6(1] 

2(1- 





6 



hi ' 






!>i- 



Following the usual prescriptions (e.g. 



Springel 2005 



Price|[2007 ), the interpolated density field is taken as 

p(x) = y^^miW{x~x.i,hi), 

i 

where tti,;, x^, h.^ are the mass, position, and smoothing 
length of the ith particle, respectively. The interpolation 
of a field variable, denoted by A, is given by 



^(x) = E 






0{h^), 



where Ai is the corresponding field property carried by 
the ith particle. Note that the density field can be seen 
as a special case of the general formula. 

Two types of pixel- wise mean for a field are calculated, 
the 



i|^d^xp(x) 



1. volume- weighted mean of the density field, 

= — V'mi / (fxW{x-:x.i,h^) 

i 

where Mi{P) — rrii Jp cfixW{:x.—Ki, hi) is the mass 
overlapping of the ith particle and the pixel, and 
the 

2. mass- weighted mean of a field (A) FM 
/pd3xA(x)p(x) 



AiP): 



/pd3xp(x) 

M{P) 



Oih'). 



To obtain a line of sight projection along the third axis, 
the pixels are chosen to extend along the third dimension, 

^^ The second line is an approximation. For the numerator, 



/ d^xA(x)p(x) = 

/■ ,3 s-^ AimiW{x-Xi,hi) 



^mjTy(x-Xj,/ij) 



resulting a two dimensional final raster image. The cal- 
culation of the overlapping Mi{P) in this circumstance 
is two dimensional. Both formulas require frequent cal- 
culation of the overlap between the kernel function and 
the pixels. An effective way to calculate the overlap is 
via a lookup table that is pre-calculated and hard coded 
in the program. Three levels of approximation are used 
in the calculation of the contribution of a particle to a 
pixel: 

1. When a particle is much smaller than a pixel, the 
particle contributes to the pixel as a whole. No 
interpolation and lookup occurs. 

2. When a particle and pixel are of similar size of (up 
to a few pixels in size) , the contribution to each of 
the pixels are calculated by interpolating between 
the overlapping areas read from a lookup table. 

3. When a particle is much larger than a pixel, the 
contribution to a pixel taken to be the center kernel 
value times the area of a pixel. 

Note that Level 1 and 3 are significantly faster than Level 
2 as they do not require interpolations. 

The rasterization of the z = 4.75 snapshot of Massive- 
Black was run on the SGI UV Blacklight supercomputer 
at the Pittsburgh Supercomputing Center. Blacklight is 
a shared memory machine equipped with a large memory 
for holding the image and a fairly large number of cores 
enabling parallelism, making it the most favorable ma- 
chine for the rasterization. The rasterization of the E5 
simulation was run on local CMU machine Warp. The 
pixel dimensions of the raster images are also listed in 
Table [T] The pixel scales have been chosen to be around 
the gravitational softening length of ^ 5 h^^Kpc in these 
simulations in order to preserve as much information in 
the image as possible. 

5. IMAGE RENDERING AND LAYER COMPOSITING 

expand, we find that 

/ d'^xWix -Xi,hi)W{yi - Xj,hj) 

= W(C,j-x,,hj) y cPxW{x-x„K) 

= [W{xi -Xj,hj) + W'{xi - Xj,hj){^ij -Xi) 
+ Omj - x,)']] 
X / d^xVy(x - Xi, hi) 

= VF(xi-Xj,/ij) / d^xWix-Xi,h,) + 0{h'^). 

Both W' and ^ij — x^ are bound by terms of 0{hj), so that the 
extra terms are all beyond 0(/i^) (the last line). Noticing that 
Pi = X], m,jW{x.i — Xj, hj) + 0{h^), the numerator 

I d^xA(x)p(x) 

= y2^i'^^ I d^xW{yi~Xi,h{) + 0{h^) 



If we apply the mean value theorem to the integral and Taylor 



The rasterized SPH images are color-mapped into 
RGBA (red, green, blue and opacity) layers. Two modes 
of color-mapping are implemented, the simple mode and 
the intensity mode. 

In the simple mode, the color of a pixel is directly 
obtained by looking up the normalized pixel value in a 
given color table. To address the large (several orders of 
magnitude) variation of the fields, the logarithm of the 
pixel value is used in place of the pixel value itself. 

In the intensity mode, the color of a pixel is determined 
in the same way as done in the simple mode. However, 
the opacity is reduced by a factor /,„ that is determined 
by the logarithm of the total mass of the SPH fluid con- 
tained within the pixel. To be more speciflc. 






log M < a 


1 

log M — a 
b-a 


log M > b 

7 

otherwise, 



where a and b are the underexposure and overexposure 
parameters: any pixel that has a mass below IC^ is com- 
pletely transparent, and any pixel that has a mass above 
10^ is completely opaque. 

The RGBA layers are stacked one on top of another to 
composite the flnal image. The compositing assumes an 
opaque black background. The formula to composite an 
opaque bottom layer B with an overlay la yer T into the 
composite layer C is (Porter fc Duff|[r984l 



C = aF + {l- a)T, 

where C, B and T stand for the RGB pixel color triplets 
of the corresponding layer and a is the opacity value of 
the pixel in the overlay layer T. For example, if the 
background is red and the overlay color is green, with 
a = 50%, the composite color is a 50%-dimmed yellow. 

Point-like (non SPH) particles are rendered differently. 
Star particles are rendered as colored points, while black 
hole particles are rendered using circles, with the radius 
proportional to the logarithm of the mass. In our ex- 
ample images, the MassiveBlack simulation visualization 
used a fast rasterizer that does not support anti-aliasing, 
whilst the fra mes of E5 are rendered using matplotlib 
( Hunter||2007 1 that does anti-aliasing. 

ihe choice of the colors in the color map has to be 
made carefully to avoid confusing different quantities. 
We choose a color gradient which spans black, red, yel- 
low and blue for the color map of the normalized gas 
density field. This color map is shown in Figure [31 Com- 
posited above the gas density field is the mass weighted 
average of the star formation rate field, shown in dark 
blue, and with completely transparency where the field 
vanishes. Additionally, we choose solid white pixels for 
the star particles. Blackholes are shown as green circles. 
In the E5 animation frames, the normalization of the gas 
density color map has been fixed so that the maximum 
and minimum values correspond to the extreme values of 
density in the last snapshot. 

6. PARALLELISM AND PERFORMANCE 

The large simulations we are interested in visualizing 
have been run on large supercomputer facilities. In or- 
der to image them with sufficient resolution to be truly 
useful, the creation of images from the raw simulation 



Figure 3. Fiducial color-map for the gas density field. The colors 
span a darkened red through yellow to blue. 

data also needs significant computing resources. In this 
section we outline our algorithms for doing this and give 
measures of performance. 

6.1. Rasterization in parallel 

We have implemented two types of parallelism, which 
we shall refer to as "tiny" and "huge", to make best use 
of shared memory architectures and distributed mem- 
ory architectures, respectively. The tiny parallelism is 
implemented with OpenMP and takes advantage of the 
case when the image can be held within the memory of 
a single computing node. The parallelism is achieved 
by distributing the particles in batches to the threads 
within one computing node. The raster pixels are then 
color-mapped in serial, as is the drawing of the point-like 
particles. The tiny mode is especially useful for interac- 
tively probing smaller simulations. 

The huge version of parallelism is implemented using 
the Message Passing Interface (MPI) libraries and is used 
when the image is larger than a single computing node 
or the computing resources within one node are insuffi- 
cient to finish the rasterization in a timely manner. The 
imaging domain is divided into horizontal stripes, each of 
which is hosted by a computing node. When the snap- 
shot is read into memory, only the particles that con- 
tribute to the pixels in a domain are scattered to the 
hosting node of the domain. Due to the growth of cos- 
mic structure as we move to lower redshifts, some of the 
stripes inevitably have many more particles than oth- 
ers, introducing load imbalance. We define the load im- 
balance penalty rj as the ratio between the maximum 
and the average of the number of particles in a stripe. 
The computing nodes with fewer particles tend to finish 
sooner than those with more. The color-mapping and 
the drawing of point-like particles are also performed in 
parallel in the huge version of parallelism. 

6.2. Performance 

The time spent in domain remapping scales linearly 
with the total number of particles N, 

Tremap ^ 0(iV). 

The time spent in color-mapping scales linearly with the 
total number of pixels P, 



T, 



color 



0{P). 



Both processes consume a very small fraction of the total 
computing cycles. 

The rasterization consumes a much larger part of the 
computing resources and it is useful to analyze it in more 
detail. If we let n be the number of pixels overlapping a 
particle, then n = K~^N~^P, where K~^ is a constant 
related to the simulation. Now we let t{n) be the time 
it takes to rasterize one particle, as a function of the 
number of pixels overlapping the particle. From the 3 
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Figure 4. MassiveBlack Simulation Rasterization Rate. We show 
the number of pixels fixed as a function of resolution (pixel scale) . 
The rate peaks at KCjT at the high resolution limit and ap- 



proaches KC2 
is not explored. 



as the resolution worsens. The n <^ 1 domain 



Table 2 

Time Taken to Rasterize the MassiveBlack Simulation . Here TV 

is the number of SPH particles, Res is the resolution (pixel scale), 

n is mean number of pixels that overlap each particle, r] is the 

load unbalance penalty averaged over patches, and the final 
column. Rate, is the number of kilopixels rasterized per second. 



Pixels Res 

(Kpc/px) 



CPUs Wall- i 
time 
(hours) 



Rate 

(K/sec) 
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58.4 
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128 


1.63 


1.45 
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22.5G 
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256 
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1.66 
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1300 


512 


3.35 


1.57 
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14.6 


1300 


512 


3.06 


1.39 


23.3 


360G 


7.3 


5300 


512 


7.65 


1.39 


37.2 


1160G 


4.2 


16000 


1344 


10.1 


1.56 


37.2 



levels of detail in the rasterization algorithm (Section 4) , 
we have 

{Ci, n<l 
C2n n ^ 1 ; 
Can, n > 1 

with C2 ^ Ci « C3 . The effective pixel filling rate R 
is defined as the total number of image pixels rasterized 
per unit time, 

R = P[t{n)N]-^ ^ nK{t(n)Y^ 
'fiKC^^, n < 1 



KCr^ 
KC; 



n~l 
n> 1 



The rasterization time to taken to create images from 
a single snapshot of the MassiveBlack simulation (at red- 
shift 4.75) at various resolutions is presented in Table [2] 
and Figure |4] 

The rasterization of the images were carried out on 
Blacklight at PSC. It is interesting to note that for the 
largest images, the disk I/O wall-time, limited by the I/O 
bandwidth of the machine, overwhelms the total comput- 
ing wall-time. The performance of the I/O subsystem 



shall be an important factor in the selection of machines 
for data visualization at this scale. 

7. IMAGE AND ANIMATION VIEWING 

Once large images or animation frames have been cre- 
ated, viewing them presents a separate problem. We use 
the GigaPan technology for this, which enables some- 
one with a web browser and Internet connection to ac- 
cess the simulation at high resolution. In this section we 
give a brief overview of the use of the GigaPan viewer 
for exploring large static images, as well as the recently 
developed GigaPan Time Machine viewer for gigapixel 
animations. 

7.1. Gigapan 

Individual gigapixel-scale images are generally too 
large to be shared in easy ways; they are too large to at- 
tach to emails, and may take minutes or longer to transfer 
in their entirety over typical Internet connections. Gi- 
gaPan addresses the problem of sharing and interactive 
viewing of large single images by streaming in real-time 
the portions of images actually needed by the viewer of 
the image, based on the viewers current area of focus 
inside the image. To support this real-time streaming, 
the image is divided up and rendered into small tiles of 
multiple resolutions. The viewer pulls in only the tiles 
needed for a given view. Many mapping programs (e.g., 
Google Maps) use the same technique. 

We have uploaded an example terapixel imagj^of the 
redshift z = 4.75 snapshot of the MassiveBlack simu- 
lation to the GigaPan website, which is run as a pub- 
licly accessible resource for sharing and viewing large 
images and movies. The dimension of the image is 
1440000 X 810000, and the finished image uncompressed 
occupies 3.58 TB of storage space. The compressed hi- 
erarchical data storage in Gigapan format is about 15% 
of the size, or 0.5 TB. There is no fundamental limits 
to size, provided the data can be stored on the disk. It 
is possible to create directly the compressed tiles of a 
GigaPan, bypassing the uncompressed image as an in- 
termediate step, and thus reducing the requirement on 
memory and disk storage. We leave this for future work. 

On the viewer side, static GigaPan works well at dif- 
ferent band widths; the interface remains responsive in- 
dependent of bandwidth, but the imagery resolves more 
slowly as the bandwidth is reduced. 250 kilobits/sec is a 
recommended bandwidth for exploring with a 1024 x 768 
window, but the system works well even when the band- 
width is lower. 

An illustration of the screen output is shown in Fig- 
ure 5] The reader is encouraged to visit the website to 
explore the image. 

7.2. GigaPan Time Machine 

In order to make animations, one starts with the ren- 
dered images of each individual snapshot in time. These 
can be gigapixel in scale or more. In our example, us- 
ing the E5 simulation (Section 2.2) we have 1367 images 
each with 0.75 gigapixels. 

One approach to showing gigapixel imagery over time 
would be to modify the single-image GigaPan viewer 



image at|http://gigapan.org/gigapans/76215/ 




Figure 5. GigaPan View of the MassiveBlack Simulation at z = 4.75. The images are screen grab from the GigaPan viewer: we have left 
the magnification bar visible. The background is an overall view of the entire snapshot. We also show zooms into region around one of the 
most massive blackholes in the simulation, as shown in the right most zoom with the largest circle. The three zoom levels are 80h"^Mpc, 
Sh-^Mpc, 800h~lKpc, from left to right. 



to animate by switching between individual GigaPan 
tile-sets. However, this approach is expensive in band- 
width and CPU, leading to sluggish updates when mov- 
ing through time. 

To solve this problem, we created a gigapixel video 
strea ming and viewing sys tem called GigaPan Time Ma- 
chine (Sargent et al. 12010), which allows the user to flu- 



idly explore gigapixei-scale videos across both space and 
time. We solve the bandwidth and CPU problems using 
an approach similar to that used for individual GigaPan 
images: we divide the gigapixel-scale video spatially into 
many smaller videos. Different video tiles contain dif- 
ferent locations of the overall video stream, at different 
levels of detail. Only the area currently being viewed 
on a client computer need be streamed to the client and 
decoded. As the user translates and zooms through dif- 
ferent areas in the video, the viewer scales and translates 
the currently streaming video, and over time the viewer 
requests from the server different video tiles which more 
closely match the current viewing area. The viewing sys- 
tem is implemented in Javascript+HTML, and takes ad- 
vantage of recent browser's ability to display and control 
videos through the new HTML5 video tag. 

The architecture of GigaPan Time Machine allows the 
content of all video tiles to be precomputed on the server; 
clients request these precomputed resources without ad- 
ditional CPU load on the server. This allows scaling to 
many simultaneously viewing clients, and allows stan- 
dard caching protocols in the browser and within the 
network itself to improve the overall efficiency of the 
system. The minimum bandwidth requirement to view 
videos without stalling depends on the size of the viewer, 
the frame rate, and the compression ratios. The individ- 
ual videos in "Evolution of the Universe" (the E5 simula- 
tion, see below for link) are currently encoded at 25 FPS 
with relatively low compression. The large video tiles re- 
quire a continuous bandwidth of 1.2 megabits/sec, and a 
burst bandwidth of 2.5 megabits/sec. 

We have uploaded an example animatiorp^ of the E5 
simulation, showing its evolution over the interval be- 
tween redshift z = 200 and z — Q with 1367 frames 
equally spaced in time by 10 Myr. Again, the reader is 
encouraged to visit the website to explore the image. 

8. CONCLUSIONS 

We have presented a framework for generating and 
viewing large images and movies of the formation of 
structure in cosmological SPH simulations. This frame- 
work has been designed specifically to tackle the prob- 
lems that occur with the largest datasets. In the genera- 
tion of images, it includes parallel rasterization (for either 
shared and distributed memory) and adaptive pixel fill- 
ing which leads to a well behaved filling rate at high res- 
olution. For viewing images, the GigaPan viewers use hi- 
erarchical caching and cloud based storage to make even 
the largest of these datasets fully explorable at high res- 
olution by anyone with an internet connection. We make 
our image making toolkit publicly available, and the Gi- 



gaPan web resources are likewise publicly accessible. 
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Figure 6. GigaPan Time Machine Animation of the E5 simulation. These are screen grabs from the GigaPan Time Machine viewer. In 
the left column we show 8 frames (out of 1367 in the full animation) which illustrate the evolution of the entire simulation volume (at time 
intervals of 2Gyears. The middle panel zooms in to show formation of the largest halo through merger event. The right panel shows some 
of the history of a smaller halo. 



